Efficient Power Spectrum Estimation for High Resolution CMB Maps 
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Estimation of the angular power spectrum of the Cosmic Microwave Background (CMB) on a 
small patch of sky is usually plagued by serious spectral leakage, specially when the map has a hard 
edge. Even on a full sky map, point source masks can alias power from large scales to small scales 
producing excess variance at high multipoles. We describe a new fast, simple and local method for 
estimation of power spectra on small patches of the sky that minimizes spectral leakage and reduces 
the variance of the spectral estimate. For example, when compared with the standard uniform 
sampling approach on a 8 degree x 8 degree patch of the sky with 2% area masked due to point 
sources, our estimator halves the errorbars at £ = 2000 and achieves a more than fourfold reduction 
in errorbars at £ = 3500. Thus, a properly analyzed experiment will have errorbars at £ = 3500 
equivalent to those of an experiment analyzed with the now standard technique with ~ 16 — 25 
times the integration time. 



I. INTRODUCTION 

Cosmic Microwave Background (CMB) is a statisti- 
cally isotropic [1] and Gaussian [2] random field. If wc 
ignore secondary effects, all of the information in high 
resolution CMB maps is encoded in the angular correla- 
tion function or equivalcntly, in the angular power spec- 
trum, Ci. The angular power spectrum is widely used 
to estimate the cosmological parameters. Accurate mea- 
surements of angular power spectrum are needed for pre- 
cise estimation of cosmological parameters. 

Over the past decade, CMB power spectrum has been 
measured over a large range of multipoles, £, by various 
groups [3-7]. And more experiments are under way to 
measure the Ci on smaller scales with high accuracy [8- 
10]. Most power spectrum analyses use uniform or noise 
weighted maps. This performs reasonably well for power 
spectra that have nearly equal power in equal logarithmic 
intervals of multipoles, i.e. £ < 1000 for the CMB. For 
smaller scales, (larger £) this method is non-optimal, as 
we show in section V. CMB power spectrum estimated 
from an incomplete sky map is the underlying full-sky 
power spectrum convolved with the power spectrum of 
the mask. This leads to coupling of modes in the esti- 
mated power spectrum. For high resolution experiments 
such as ACT and SPT which will map the small scale 
anisotropies of the CMB on small patches of the sky, this 
mode-mode coupling will be a serious problem. The rea- 
son is that CMB power spectrum is very red on those 
scales (it falls off as £~ A at large £) and hence is highly 
vulnerable to the leakage of power due to mode-mode 
coupling. There are two methods to remedy this: to 
taper the map near the sharp edges, and to pre-whitcn 
the CMB power spectrum. In order to minimize the loss 
of information due to applying a taper to the map, we 
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use the multitaper method [11]. This method involves 
weighting the map with a set of orthonormal functions 
which arc space limited but maximally concentrated in 
the frequency domain. Power spectrum of each of these 
tapered maps is a measurement of the power spectrum of 
that map with a different amount of mode coupling. Fi- 
nal power spectrum is obtained from a particular linear 
combination of these tapered power spectra that mini- 
mizes the bias in the estimated power spectrum. The 
use of multiple tapers also reduces the error-bars in the 
measured power spectrum. 

Mode coupling is less harmful if the map has a nearly 
white power spectrum. Traditionally, an inverse covari- 
ance matrix weighting is used in analysis to prewhiten 
the maps [12]. This method works well, but is a com- 
putationally expensive, non-local operation and may be 
complicated to implement, specially for high resolution 
experiments [13, 14]. Wc propose a simple and local 
prcwhitening operator in real space (§ IV) that is fast 
to implement and reduces the bias due to the leakage 
of power. This method prevents the unnecessarily large 
error bars at I > 1500 due to the point source masks. 
Usually masks have sharp edges and holes at the po- 
sitions of point sources. This leads to a mode-coupled 
power spectrum that is highly biased at large I. Decon- 
volution of the mode-coupled power spectrum is a well- 
studied problem in the CMB data analysis literature [15] 
and has been applied to many experiments. But decon- 
volution of a highly biased power spectrum leads to large 
error bars in the final power spectrum at large I. The 
mode coupling problem will be worse for the upcoming 
set of CMB experiments as bright point sources will be 
much more of a limiting foreground at high resolution. 

As we show in § V, prewhitcning followed by the multi- 
taper method for power spectrum estimation reduces the 
error-bars (specially at large t) in the decoupled power 
spectrum (cf. Figs. 14 & 15). 

We begin with a review of the multitaper method in 
one-dimension in § II and discuss the salient features of 
the method, generalizing it to the two-dimensional case. 
Next, we discuss the statistical properties of multitaper 
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spectrum estimators. As a simple application, we demon- 
strate the method in context of CMB power spectrum 
estimation in § III. Next, we formulate the prewhitcning 
method (§ IV) and apply it to the case of CMB power 
spectrum estimation in presence of masks. In § V, we 
describe the algorithm for deconvolving the power spec- 
trum and the implications of the multitaper method and 
prewhitcning in its context. Wc summarize and conclude 
in 6 VI. 



A. Notations 

Throughout this paper, we will refer to spatial coor- 
dinates as the x space (this may be an angular coordi- 
nate in radians on the sky, or a comoving distance in 
/i _1 Mpc, etc.) and the reciprocal space as the k space 
(which would be the multipole space £, or the Fourier 
modes in /i _1 Mpc, etc). The continuous Fourier Trans- 
form conventions adopted here are, 



II. A BRIEF REVIEW OF THE MULTITAPER 
METHOD 



The problem of estimating the power spectrum of a 
stationary, ergodic process, sampled at discrete inter- 
vals and observed over a finite segment of its duration 
of occurrence, is an old and well-studied one (for an ex- 
tensive treatise, see [11]). Several methods have been 
traditionally used for power spectrum estimation in one- 
dimension. These include non-parametric methods like 
the periodogram, the lag-window estimators, Welch's 
overlapped segment averaging [16] and the Multitaper 
method [17], and parametric methods like the maximum 
likelihood estimation. In this paper, we generalize the 
one-dimensional Multitaper method to two-dimensions 
and adapt it to handle real data with noise and masks 
on a two-dimensional flat Euclidean patch. We discuss 
its applications specifically in the case of CMB power 
spectrum estimation. 

The most basic spectral estimation method is to take 
the square of the Fourier Transform (FT) of the observed 
data. Taking the FT of a finite segment of data is equiv- 
alent to convolving the underlying power spectrum with 
the power spectrum of a top-hat function. As the latter 
has substantial sidelobe power, it leads to spectra leak- 
age and the resulting spectrum is highly biased. Most of 
the non-parametric methods for power spectrum estima- 
tion utilize some kind of a data taper (a smooth function 
that goes smoothly to zero at the edges of the observed 
segment) to minimize the effect of spectral leakage. Such 
smoothing reduces the bias in the estimator at the cost 
of lower spectral resolution. As the taper down-weights 
a fraction of the data, one is left with an effectively lower 
sample size. Since tapering also smooths in frequency 
space, it essentially leads to a loss of information which 
is reflected in the increased variance of the final estimate. 
The first attempt at ameliorating these disadvantages of 
using a data taper was addressed in a seminal paper by 
[17] which laid down the basis of the Multitaper Method 
(MTM). The basic idea of MTM is to apply multiple or- 
thogonal tapers with optimal spectral concentration to 
minimize the loss of information due to tapering. 



F(k) = I (fx F(x) exp(-ik • x) (forward) (1) 



F(x) 
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F(k) cxp(ik • x) (inverse), (2) 



where n is the dimensionality of the space. 

For a stationary process, -F(x) the power spectrum is 
defined as, 



(27r)"P(fc) = (^(k)F(k)). (3) 



B. 1-D Multitaper Theory 

Although power spectrum estimation for the CMB is 
an inherently two-dimensional problem, we will begin by 
discussing the multitaper theory in one dimension. This 
is because the essential features of the theory are eas- 
ier to understand in one dimension and can be trivially 
generalized to higher dimensions. 

We consider a stationary, stochastic, zero mean pro- 
cess F(x) sampled at N discrete points, Xj sampled at 
regular intervals of size Ax. Let P{k) be the true under- 
lying power spectrum of the process. Our problem is to 
estimate P(k) using the sample of size N. 

The Nyquist frequency for the problem is given by 
}n = fcjv/(27r) = l/(2Ax) and the fundamental fre- 
quency by fo = fco/(27r) = l/(iVAa;). In the following, 
we will assume Aa: = 1 for simplicity. 

Let us contemplate windowing our data by some func- 
tion G(x), generating the product, 

y( Xj ) = G( Xj )F( Xj ) (4) 

and taking the power spectrum of the windowed data as, 



(27r)P(fc) = f(k)y(k) = 



N 



j'=i 



kx 



, (5) 



where y is the Fourier transform of y{x). The quantity 
P can be thought of as an estimator of P(k), such that 
its ensemble average is related to P(k) through 



P(k) 



— hp 



dk' 

W) 



T(k - k')P(k'). 



(6) 
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This means that P is an estimator of the true power 
spectrum convolved with a spectral window junction, 



T{k) = G(k) . 



(7) 



If r could be designed such that T(fc) = (2iv)8(k) then P 
would be an exact unbiased estimator of P. However, a 
function like G which is spatially limited in extent cannot 
be arbitrarily concentrated in the frequency space. If the 
window function is a top-hat, its power spectrum will 
be a sine 2 function[32] with substantial sidelobes. This 
will lead to the aliasing of power on various scales, an 
effect known as spectral leakage or mode coupling. Mode 
coupling is specially damaging for a spectrum which is 
highly colored or structured. 

The multitapcr method (MTM) consists of finding 
a set of orthogonal window functions or tapers, which 
are maximally concentrated in some predetermined fre- 
quency interval. With the set of tapers, one can gener- 
ate several approximately uncorrelated estimates of the 
power spectrum. This is superior to the plain Fourier 
Transform (Pcriodogram) because it not only attempts 
at remedying mode coupling errors but also helps de- 
crease the uncertainty in the estimated power spectrum 
by generating independent realizations of the same power 
spectrum with information from different section of the 
data. We formulate the method below. 

We desire a set of tapers, such that each of them is 
spatially limited, Gj (j = l,...,N) and has its power 
r(/c) optimally concentrated in some frequency interval, 
k G [-2nW 7 2irW], with W < f N = k N /(2n). Here we 
have introduced the shorthand notation Gj for G{xj). 
Concentration is quantified by the following quantity, 



(3 2 (W) 



I-MV T ( k ) dk 



(8) 



which is basically the fractional power of the taper inside 
the desired interval. Remembering that, 



N 



(9) 



the above equation can be re-written as, 



j'=lj=l / 3=1 

(10) 

It is easy to see that the sequence Gj that will maximize 
(3(W) must satisfy, 



h x <i - r) 



G y = \ a (N,W)Gj, (11) 



for j = 1,...,N. This can be immediately recognized as 
an eigenvalue problem, 



where A is the N x N Tocplitz matrix, 
_ sm[2nW(j ~ j')] 



(13) 



and G is the vector, index limited from j = 1 to j = N 
which has the highest concentration in the frequency 
interval [— W, W]. Here a denotes the indices of the 
different eigenvalues of the problem. The solution to 
this eigenvalue problem is well known [18]. There are 
TV nonzero eigenvalues of the problem denoted by X a 
(a = 0,1, (JV — 1)) with corresponding eigenvectors 
v Q . The elements of each of the N eigenvectors consist 
of a finite subset of the discrete prolate spheroidal se- 
quence (DPSS). The zeroth eigenvector v°(iV, W) which 
has the highest eigenvalue Ao is composed of the zeroth 
order DPSS, the eigenvector v 1 (A r , W) having eigenvalue 
Ai < Ao is composed of first order DPSS sequence and so 
on. 

Some salient properties of the N eigenvectors and 
eigenvalues are as follows, 

1. The eigenvalues are bounded by and 1: 

< v a < 1. 

2. The eigenvectors are orthogonal and can be stan- 
dardized so that they are orthonormal, 
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3. The eigenvectors form a basis for an iV-dimensional 
Euclidean space. 

4. Usually the eigenvectors are ordered according to 
decreasing eigenvalues. The first 2^1^—1 eigenval- 
ues are close to unity (most concentrated) and the 
eigenvalues rapidly fall to zero thereafter. This be- 
havior is illustrated in Fig. 1. The number 2NW—1 
is often referred to in the MTM literature as the 
Shanon Number. 

Onc-dimcnsional DPSS taper generation algorithms are 
usually included in standard signal processing softwares 
(e.g. the dpss module of Matlab). For the purpose of 
this paper, we used a Fortran 90 implementation of the 
original algorithm by [19]. Examples of DPSS tapers 
and their corresponding spectral window functions are 
displayed in Fig 2, where the gradual worsening of the 
leakage properties of the tapers are apparent. 
The Fourier transforms of the tapers, 



(14) 



also have the interesting properties: 



1. They are orthonormal over the frequency range 



AG = X a (N, W)G 



(12) 
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FIG. 1: Eigenvalues corresponding the different orders of 
DPSS tapers. Two cases with NW = 3 and 6 are shown for 
N — 50. Spectral concentration of the tapers rapidly worsen 
beyond a = 2NW - 1. 



2. They are also orthogonal over the frequency do- 
main -2itW < k < 2nW, 



2kW , h 

—v a *(k)v^(k) = X a S aP . (15) 



2irW 



2ir 



This means that the functions v q /VAq form an 
orthonormal set on the inner interval —2nW < k < 
2ttW. 

Having generated the eigentapers, we can form M ap- 
proximately uncorrelated estimators of the power spec- 
trum with the first M eigenvectors having the best con- 
centration, 



N 
3=1 



(16) 



where a = 0, M — 1. 

Then we can form a weighted mean of the tapered 
power spectra, often called the eigenspectra, to generate 
the final estimate of P(k), the simplest form of which is, 



P 



MTM 



E 



M-l 



A, 



(17) 



It can be shown that this is the optimal way of estimating 
P in the case where the process is white noise. For col- 
ored spectra, a more sophisticated approach is required, 
which leads to the adaptive multitaper method (AMTM) 
to be discussed in the following section. 

The important point to note here is that the MTM 
or variants of it, aim to restore the information lost to 
a single taper algorithm by weighting different parts of 
the data by an orthogonal set of tapers, thereby reducing 
the variance in the final estimate. Lower variance comes 
at the cost of decreased spectral resolution. One should 



bear in mind that the choice of the resolution W is com- 
pletely dependent on the analyst. Remembering that the 
number of useful tapers is (2NW — 1) and a better spec- 
tral resolution, i.e. smaller W means that there will be 
fewer tapers to work with. The choice of W will in most 
cases be dictated by the type of and the features in the 
power spectrum being estimated. 



C. Adaptive MTM 

As noted above, if the underlying spectrum is white 
and only modest spectral resolution is needed in the anal- 
ysis, many eigenspectra can be simply combined with 
scalar weights to get a good estimate of the true spec- 
trum. However, this is not the case for spectra which are 
colored and have large dynamic range. For example, the 
Cosmic Microwave Background (CMB) power spectrum 
C t falls like £" 4 beyond I ~ 1000. In cases such as this, 
only the first few tapers are good at avoiding aliasing of 
power due to mode coupling. As more and more tapers 
are used, the estimated power spectrum gets more and 
more biased. 

The adaptive multitaper method (AMTM) aims at 
mitigating this problem, thereby allowing the use of a 
larger number of tapers even for a colored spectrum. In 
the following we briefly sketch the AMTM method. 

According to the Cramer spectral representation of a 
stationary process [20], a stationary zero mean process 
can be represented as, 



F[x) = / e lkx dZ{k) 
-fe» 



(18) 



for all x, where dZ is an orthogonal incremental process 
[21, 22]. The random orthogonal measure dZ(k) has the 
properties, 



(dZ(k)) = 0; ( \dZ(k)\ 2 ) = P{k)dk 



(19) 



where P(k), as before, is the true underlying spectrum. 

The Fourier transform of the data weighted by an 
eigentaper can be written as, 



A' 



y a (k) 



\Xi e 



^ — ik : 



Xycain- 

d a (k- k')dZ(k') 



(20) 



(21) 



using (18). Note that this contains information from the 
entire Nyquist range. 

Now consider the case where the signal F(x) is con- 
volved with a perfect bandpass filter from k — 2irW to 
k + 2irW to yield the unobservable yet perfect eigencom- 
poncnt, 



y a (k) 



k+2-ixW 



i a (k-k') 



dZ(k') (22) 



k'=k-2-KW 



5 





0.4 














0.2 
0.0 
-0.2 
-0.4 












~> 


Order = 


Order = 2 


Order = 4 


Order = 6 


Order = 8 



10 20 30 40 10 20 30 40 10 20 30 40 10 20 30 40 10 20 30 40 

Real space cooridinate (x) 




-0.4-0.3 0.0 0.2 0.4 -0.4-0.2 0.0 0.2 0.4 -0.4-0.2 0.0 0.2 0.4 -0.4-0.2 0.0 0.2 0.4 -0.4-0.2 0.0 0.2 0.4 

Fourier space frequency (f) 



FIG. 2: Examples of DPSS tapers and the corresponding spectral window functions for the case N = 50 and NW = 6. Upper 
panel: Real space form of the tapers of orders 0, 4, 6 and 8. Lower panel: The spectral window functions corresponding to 
the tapers in the upper panel. For simplicity, the window functions are shown only in the range — 3W < f < 3W of frequency. 
The vertical dotted lines denote the edges of the bandwidth ( — W, W) within which the tapers are designed to be optimally 
concentrated. Tapers are ordered such that spectral leakage progressively increases for tapers of higher order. 



which contains information only from the interval 
(— 2ttW, 2itW). Note that in accordance with (15) we 
have used the correct orthonormal form for FT of the 
tapers inside this interval. 

The quantity y(k), although fictitious in the sense that 
it cannot be computed from the data, has the desirable 
property that, 



P{k')dk' (23) 



i.e. it is the unbiased estimator of the true power spec- 
trum smoothed by a strict bandpass filter of width AirW. 
Therefore, in a multitaper setting, when combining dif- 
ferent tapers with weights, the weights should be cho- 
sen such that the eigenspectrum obtained by each of the 
weighted tapers is as close as possible to this ideal esti- 
mate. This forms the basis of the AMTM, where one re- 
places the scalar weights by frequency dependent weight 
functions b a {k) which minimize the mean squared error, 



y(k) 


9 \ r k+2-!rW 

)- 


'\v a (k-k')\~ 




y/X a 




1 Jk'=k~1-nW 



MSE«(fc) = ( y a (k)~b a (k)y a (k) 



(24) 



The minimization leads to the following expression for 



b a (k) = 



KP{k) 



X a P(k) + (1 - X a )a 2 



(25) 



where, a 2 = J ^ ||P(fc) = varfx^}, the variance of the 
map, using Parseval's theorem. 



With this the minimum mean square error becomes, 
P(k)(l-X a )a 2 



MSE Q (fc) 



X a P(k) + (1 - X a )a 2 



(26) 



and the best estimate of the power spectrum has the 
form, 

pAMTM (fc) = bl(k)P a (k) 
E«=obl(k) 

where P a is the single taper power spectrum estimate 
with the a th taper i.e. the eigenspectrum of order a. In 
case of a white noise spectrum P(k) = a 2 , and therefore 
b a {k) = y/X^, which gives us the formula (17) for the 
simple MTM. Also, in this case MSE Q = (1 - X a )a 2 . 
This shows again that for white noise, for which the first 
2NW— 1 tapers have eigenvalues close to unity, the mean 
squared error is negligible. 

Note that the estimation of the power spectrum re- 
quires the evaluation of the optimal value for the weights 
b a (k), which assumes knowledge of the true power spec- 
trum. But the latter is precisely what we are trying to 
estimate. Therefore, this method has to be implemented 
iteratively according to the following steps, 

1. We use the first one or two tapers (having the least 
spectral leakage) to form a first estimate of the 
power spectrum P(k) via (17). 

2. We use (25) to estimate the weights b a (k). 

3. With the b a 's estimated, we use (27) with M < 
2NW — 1 tapers to get the second estimate of the 
power spectrum. 
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4. Using this new estimate, re-evaluate the weights 
b a (k) and repeat steps 2-4. 

In the following, we will generalize the AMTM to two 
dimensions. 



D. AMTM in two dimensions 

So far we have only considered onc-dimcnsional tapers. 
In the signal processing literature, generalizations of the 
multitaper method to higher dimensions have not been 
widely discussed. Some of the early works include [23- 
25]. Relatively recently, straightforward generalizations 
of the method to higher dimensional cuclidcan spaces [26] 
and to patches on the surface of the sphere [27-29] have 
been formulated. In view of the upcoming CMB experi- 
ments (that will map the sky on small scales at high res- 
olution), we will be mostly interested in the power spec- 
trum estimation on regularly sampled two-dimensional 
flat spaces like a projection of a small patch on the sky. 
As such we will discuss the 2-D extension of the multita- 
per method as discussed in [26]. 

Two-dimensional tapers are constructed from outer 
products of one-dimensional tapers. We assume a two- 
dimensional map spanned by the coordinates (x 1 ,^; 2 ) . 
The data should be sampled regularly, but the sampling 
intervals (pixel sizes), Ax 1 , and the number of pixels in 
each direction, Ni, can be different (i = 1,2). Let us re- 
label the one-dimensional taper of the previous section as 
v Q;JV where the extra superscript denotes the number of 
pixels in a given direction. Then a two-dimensional taper 
of order (a±; a 2 ) and size JVi x N 2 can be constructed out 
of the outer product of two one-dimensional tapers: 



N X N 2 ) 



r {a 2 ,N 2 ) 



(28) 



where we have treated v column vector. The spec- 
tral concentration eigenvalue of a two-dimensional taper, 
is easily shown to be the product of the eigenvalues of 
the one-dimensional tapers out of which the former is 
constructed: 



\a 1 a 2 ;N 1 N 2 ) ~ ^a 1 M 1 ^a 2 ,N 2 - 



(29) 



Two-dimensional tapers constructed in this way inherit 
most of the properties of one-dimensional tapers such as 
orthogonality on the sample plane, and optimal spectral 
concentration and orthonormality in the frequency plane. 
The two-dimensional power spectrum estimator corre- 
sponding to these tapers is defined, similar to the one- 
dimensional case, as a weighted sum of approximately 
independent tapered power spectra 



P AMTM (k 1 ,k 2 ) 



E ai ,q 2 6 2 aia2) (fc 1 ,fc 2 )P^"^(fc 1 ,fc 2 ) 

V b 2 : 



(30) 

where we have dropped the NiN 2 portion of the labels 
for simplicity. Note that in the above formula the weights 



b( ai a 2 ) depend on the eigenvalues \r aia2 ) and are given 
by an equation analogous to (25) and are to be esti- 
mated itcratively. The quantity (2tt) 2 P^ aia ^ (k ±1 k 2 ) is 
the eigenspectrum of order (aia 2 ) and is given by 



$(aia 2 ) _ 



FT 



(31) 



where M = {M(xi, Xj)} is the two-dimensional map, the 
power spectrum of which is being estimated. 

We will introduce a bit of notation at this point. We 
will designate the two parameters that control the mul- 
titaper method as: 

• Ntap- The number of tapers used. Its value will 
be written in the form M 2 , where M denote the 
number of one-dimensional tapers. For example if 
2 one-dimensional tapers are used to create 4 two- 
dimensional tapers via outer products, then we will 
quote N tap = 2 2 . 

• N res : We will use this as the shorthand for the res- 
olution parameter NW. For example, N res = 3 
will mean that the half-bandpass chosen for gen- 
crating the tapers is three times the fundamental 
frequency. 



E. Statistical Properties 

Statistical distribution of the power spectra of a Gaus- 
sian random field realized on a map with periodic bound- 
ary conditions, such as a full sky CMB field, can be de- 
scribed in terms of the simple analytic distributions. This 
property stems from the fact that each Fourier-mode (or 
spherical harmonic component) of such a map is a stati- 
cally independent quantity. On a finite patch of the sky 
or for any map with a non-periodic boundary condition, 
these modes get entangled due to the convolution with 
the window, and are no longer amenable to such sim- 
ple descriptions. As we will discuss in this subsection, 
the multitaper method approximately restores many of 
these nice properties of random fields on a finite patch, 
and makes the statistical properties of the spectral esti- 
mators describable via simple and intuitive analytic ex- 
pressions. 

Most of the statistical properties of the different spec- 
tral estimators that we have discussed so far stem from 
the basic result that for most stationary processes with 
a power spectrum P(k) that is continuous over the in- 
terval [— fcjvi fcjv], the simple FFT power spectrum (pe- 
riodogram) P PM (k) is distributed as, 



P FM {k) = 



d \P{k)xl/2, for < k < k N ; 



\P(k)xl 



for k — or k = k 



N- 



(32) 



asymptotically as N — > 00. Here = means "equal in dis- 
tribution" , which means that the statement 11 X = ax 2 " is 
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FIG. 3: Top panel: The first 4 two-dimensional tapers on a NxN grid with N — 50 and NW = 5. The order of the one- 
dimensional tapers corresponding to each taper is indicated on the top. Bottom panel: The logarithm of the spectral window 
functions F(fc 1 , k 2 ), corresponding to the tapers on the upper panel. The color scale on the lower panel are standardized to be 
(1CT 30 , 1) times the maximum in each plot. The white dotted lines represent the location of ±2nW wavenumbers. 



equivalent to saying that the random variable X has the 
same distribution as a chi-square random variable with 
v dcgrccs-of-frecdom (dof) that has been multiplied by 
a constant a. For a Gaussian white noise process the 
above result is exact for any N. Also, for the asymp- 
totic case, the power spectra for two different frequencies 
k and k' are uncorrelated. Although these results are 
true for asymptotically large N, in a finite N case, they 
approximately hold for the N/2 + 1 independent Fourier 
frequencies kj = (2ivj)/(NAx), if N is large enough. 

Using the above result, it is possible to predict the ap- 
proximate distribution of an AMTM spectrum estimator 
(30). Using an equivalent dcgrees-of- freedom argument 
(see appendix A for details), it can be shown that for 
N -> oo, 



P 



AMTM 



(k) 



d(P AMTM (k)) 



,(k) 



where, 



i/(k) = 



2 (j2 ai ,a 2 k(aia 2 ) 0*) 

V hf fkT 

Z— *ai,a2 (ckiQ2)^ ' 



(33) 



(34) 



Therefore, for finite N, it is reasonable to expect that the 
AMTM power spectrum at each pixel is approximately 
distributed as ^p AMTM ^ y/^/v with v given by the above 
equation. Note that in the case of MTM, b 2 aiao = \ ai a 2 - 
If we are only using tapers with eigenvalues close to unity, 
then i/(k) ~ 2M, M being the total number of tapers 
used. 

Now we turn to the approximate form of the distribu- 
tion for the power spectrum after it is binned in annular 



rings in k space, which we denote by P B {k b ), 

P B (h) = ^- b J2P AU ™(h,k j ) (35) 



i,jeb 



where the sum is over all pixels that fall inside bin b 
and Nb is the number of observations in that bin. Us- 
ing a similar argument as before, it can be shown (see 
appendix A) that, 



P B (h) 



where the degree-of-freedom, v b is given by, 
1 2JV? 



res \ ^ 

,.(L. : 



'"• N b 2 .^Kfc^r 



(36) 



(37) 



with i/(k) given by (34). If the degree-of-freedom vari- 
able is also slowly varying, then this implies v b = 
Nf,/(2N 2 es ) f(|k| ~ kb), which is the expected result 
for the sum of Nb/(2N^ es ) independent identically dis- 
tributed xt variables. Note that the appearance of the 
iV^es factor essentially arises from the fact that AMTM 
significantly correlates N^ es nearby pixels, and therefore 
only Nb/(2N^ es ) "super-pixels" are approximately statis- 
tically independent. 

In case of the periodogram, v = 2 and N res = 1, so 
that Vb — Nb- In case of MTM with M tapers with 
good leakage properties, v ~ 2M and therefore, Vb — 
N h M/N* 
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FIG. 4: One realization of the CMB map on a 192 x 192 pixel 
grid. The physical size is 8 degrees on a side. Estimation 
of the power spectrum is done on the 4 degree x 4 degree 
subarea indicated by the rectangle. The color scale represents 
temperature fluctuations in micro-Kelvin. 



III. APPLICATION TO THE IDEAL CMB MAP 

In this section, we will illustrate the multitapcr method 
by applying it to flat-sky cosmic microwave background 
(CMB) maps. We shall assume the maps to be pure CMB 
without noise or masks and with uniform weight. We also 
assume a square map with square pixels for simplicity, 
although all the results shown hold for non-square maps 
and/or non-square pixels. 

We adopt a deliberate change of notation at this point 
to make the following sections more compatible with ex- 
isting CMB literature. Since we will dealing with angular 
co-ordinates on the sky, we call the real space variable 9 
rather than x. We also call the dual Fourier space, the 
£ space, instead of the k space. The latter is motivated 
by the fact that the full sky CMB power spectra are ex- 
pressed in terms of multipole moments, which are de- 
noted by £, and the k vector is the correct generalization 
of the £ modes in the flat-sky case. We will also denote 
the power spectrum P{k) by Cg. 

We generate the CMB maps as Gaussian random real- 
izations from a theoretical power spectrum Cg. In order 
to simulate the effect of non-periodic boundary condi- 
tions we first generate larger maps from which the de- 
sired region is extracted. In the present example, we 
generated 5000 realizations of 192 x 192 pixel CMB map 
having a physical size of 8 degrees on a side. This im- 
plies that the sampling intervals (i.e. pixel scales) along 
cither axes is 59 — 2.5'. We extract a 96 x 96 pixel 4 



degree square sub-map from the center of each such re- 
alization and perform the power spectral estimation on 
this sub-area (see Fig. 4). Given the physical size and 
the number of pixels, one finds the Nyquist and funda- 
mental values of £ to be £^ vq = ir/(A8) = 4320 and 
£fund = (2tt)/4° = 90. Therefore, if we choose a resolu- 
tion parameter of N res = 2 the half-bandwidth outside 
which spectral leakage is minimized is We = 180. 

We perform the straight FFT or periodogram (labeled 
PM), MTM and AMTM power spectrum estimation on 
5000 random realizations of the CMB map. The results 
are shown in Fig 5. These plots show the mean power 
spectrum binned uniformly in £ in bins of A£ = 180. The 
error bars correspond to the 2a spread in the distribu- 
tion of their values. Two features are immediately ap- 
parent from this figure. First, it shows that the PM has 
significant spectral leakage and produces a power spec- 
trum which is highly biased. Although the eigenvalue 
weighted MTM is better in this respect, it still suffers 
from significant bias at high £'s because the higher or- 
der tapers cannot guard efficiently against spectral leak- 
age. The AMTM seems to do very well in minimizing 
spectral leakage and producing an approximately unbi- 
ased estimate of the input power spectrum. Second, the 
2a spread in the uncertainty of the binned value in the 
PM case is much higher than the corresponding spread 
for the AMTM. This is one of the main reasons for per- 
forming tapered power spectrum estimation with multi- 
ple tapers, as has been stressed before. As discussed in 
§ V this property becomes extremely important when the 
tapered power spectra are decorrelated via a MASTER- 
like [15] algorithm. The nearly unbiased nature of the 
AMTM power spectrum translates to errorbars in the 
deconvolved spectrum which are several factors smaller 
than those obtained from a periodogram in the large £ 
regime. 

Now, we will compare the distributions of the various 
power spectrum estimators with the theoretical expecta- 
tions of section II E. For this purpose, we choose three 
multipoles £ ~ 1000, 2000, 3000 at which we study both 
the single pixel and the binned probability distributions. 
The reason for choosing these three numbers is that the 
first multipole is an example where all three methods are 
approximately unbiased, the second multipole is a case 
where the PM is biased but the MTM and AMTM are 
almost unbiased, while the third multipole represents a 
regime where only the AMTM is close to an unbiased 
estimate. 

For the single pixel case, we choose three pixels P) 
on the -£-space map, such that values of \£\ are close to 
the three multipoles as described above. We store the 
values of Ci realized in the 5000 Monte Carlo simula- 
tions in these pixels, and draw up the probability distri- 
bution of the quantity £{£ + 1)C i / {2tt) from these values, 
for each of the three methods. The results are shown 
in Fig 6. In the top panel, we show the results of the 
same experiment that was used to generate Fig. 5, i.e. 



with N r 



2.0 and N, 



tap 



3 . For the leftmost plot, 




FIG. 5: Comparison of the simple periodogram method (labeled PM), the eigenvalue weighted multitaper (labeled MTM) 
and the adaptive multitaper (labeled AMTM) methods for estimating the power spectrum of a CMB map. In each plot, the 
continuous line represents the theory power spectrum used as an input for the Monte Carlo simulations. The open circles 
represent the mean values in each I bin, averaged over 5000 realizations, while the vertical lines represent the 2a spread. The 
bin width for this figure was A£ = 180. For the MTM and AMTM methods N tap = 3 2 tapers with N res = 2 were used (see text 
for details). Note that the power spectra for the multitaper methods appear smoothed because they are convolved with the 
window function of an effective taper. Standard decorrelation techniques, like the MASTER algorithm [15] can be employed 
to de-bias and deconvolve all the above power spectra, but in the first two cases, where the mode-coupled power spectra are 
biased, decorrelation leads to bigger uncertainties in the deconvolved power spectra (see § V). 



all three methods of power spectrum estimation are ap- 
proximately unbiased. The MTM and AMTM methods 
are almost identical here, because the spectrum does not 
have a large local slope at this multipole so that even in 
the adaptive method, the higher order tapers do not need 
to be down-weighted to reduce spectral leakage. There- 
fore, the equivalent degrees-of-freedom come out to be 
the same and ~ 2iV tap = 2 x 3 2 = 18, as expected for 
a white spectrum from (34). For the middle plot, the 
multipole t is such that the local spectrum is moder- 
ately colored and spectral leakage for a periodogram be- 
comes apparent in the biased estimate of the mean power 
spectrum it generates, as indicated by the dashed verti- 
cal line. The AMTM remains approximately unbiased 
while the MTM is only slightly biased, but the equiv- 
alent degrees-of-frcedom for AMTM is now lower than 
that of the MTM method. This is expected because to 
reduce spectral leakage the weights associated with the 
higher order tapers have been reduced in the AMTM 
spectrum, thereby lowering the degrees-of-freedom. This 
trend continues to higher multipolcs, and as shown in 
the rightmost plot of in the figure, the spectral leakage 
is kept at a minimum only in the AMTM power spec- 
trum by heavily down- weighting the higher order tapers, 
while both MTM and PM become highly biased. An im- 
portant point to note here is that not only does AMTM 
minimize the spectral leakage (and hence bias) but it also 
takes advantage of multiple tapers by reducing the error 
in the estimate. This is most easily seen by comparing 
the width of the distributions of AMTM spectra against 
that of the PM spectra. The lower panel of Fig. 6 es- 



sentially shows the same features but for N res = 3.0 and 
Ntap — 5 2 . A comparison of the upper and lower panels 
also illustrate how by using more tapers the uncertainty 
in the power spectrum can be reduced. We would like 
to remind the reader once again here that the number of 
usable tapers depend on the spectral resolution chosen - 
poorer spectral resolution (in this example N res = 3 vs. 
2) allows the use of more tapers (N tap = 5 vs. 3) and 
consequently an estimate with lower error. 

Now, we turn to the distribution of the bandpowers 
C , i.e. the mean power spectrum in annular bands in 
the t space. This is illustrated in Fig. 7. The top two 
panels correspond to the same multitaper parameters as 
in Fig. 6, binned uniformly in t space with bins of width 
A£ = 180, while in the third panel we show the case for 
another experiment with N res = 4.0 and Ntap = 4 binned 
at A£ = 360. The curves plotted over the points are 
the approximate theoretical distributions expected from 
(37). We note that the theoretical curves are always a 
good fit for the AMTM method. This is mainly because 
of the fact that the AMTM is approximately unbiased 
at all multipolcs, while the MTM or PM become biased 
except at the lowest multipolcs. Since near-unbiascdness 
was an assumption adopted in deriving the binned dis- 
tributions, the latter two methods suffer from mismatch 
with theory wherever they are highly biased. For the 
N res = 4.0 case, we chose a bin width of 360 so as to 
make it large enough to avoid splitting a super-pixel be- 
tween bins. Also, note that the distribution of the binned 
power spectra are close to Gaussian, an expected result, 
as binning essentially involves combining many variables 
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with almost identical distributions. 

Next, we investigate the covariance between different 
bins, in order to study how the bandpowers are correlated 
with each other. We define the scaled covariance matrix, 

Cab = , ab (38) 

V^aa^bb 

where, 

C ab = ((C a - (Ca)) 2 (C b - (C b )) 2 ) . (39) 

We use our Monte Carlo simulations to estimate the 
quantity above. The results arc displayed in Fig. 8. 
These figures illustrate an extremely desirable feature 
of the AMTM estimator i.e. the bandpowers are essen- 
tially uncorrelated beyond the spectral resolution A£w = 
2 N res ifund set by the taper parameters. This implies 
that if we set the bin widths to be same as A£w then ad- 
jacent bins will be uncorrelated. If our bins arc smaller, 
they will be correlated through a mode coupling matrix 
which is fairly diagonal, and hence can be easily inverted 
to decouple them, an issue we will briefly touch upon in 
§ V. 

In the following section, we turn to the more practical 
issue of dealing with CMB maps with noise and point 
source masks in them. 



1. Perform local real space convolution of the map 
with designed kernels so as to make its power spec- 
trum as flat (white) as possible, at least over the 
range of multipoles which is affected most by alias- 
ing of power due to the point source mask. We refer 
to this procedure as "Prewhitcning" . 

2. Perform an AMTM power spectrum estimation of 
this prewhitcned map in order to minimize any ad- 
ditional aliasing of power due to sharp edges of the 
map. 

3. As the prewhitcning operation was a convolution 
whose Fourier space form is (preferably analyti- 
cally) known, divide the power spectrum of the 
prewhitcned map with the Fourier space form of the 
prewhitening function to recover the power spec- 
trum. 

Note that the design of the prewhitcning function will 
be specific to the type of signal being considered. In the 
following, we will demonstrate the prewhitcning of CMB 
maps, with forms of the prewhitcning function specific to 
the features in the CMB power spectrum. We will first 
consider a noiseless map, in order to motivate the basic 
form of the prewhitening operation and then generalize 
to maps with white noise. 



IV. PREWHITENING FOR CMB MAPS WITH 
MASK AND NOISE 

The maps of the cosmic microwave background pro- 
duced by any experiment will invariably contain instru- 
mental noise and regions, like bright point sources, that 
are usually masked out before estimating the power spec- 
trum. The raw sky map also contains other astrophysical 
contaminants which we will neglect for the purpose of this 
discussion. 

Application of point source mask to a map is essen- 
tially replacing the pixel values in an area containing each 
point source with zeros. Masking is therefore equivalent 
to multiplying the map with a function which is unity 
everywhere except inside discs of varying sizes, where 
it is zero. This can also be thought of as successively 
multiplying the map with a mask for each point source. 
Taking power spectrum of the final masked map is there- 
fore equivalent to successively convolving the true power 
spectrum with power spectra of a series of such singlc- 
point-source masks. As the size of the discs will in general 
vary for each such function, the true power spectrum will 
be convolved with functions that have power over various 
ranges of multipoles. Although the multiplication with 
tapers will guard against aliasing of power due to sharp 
edge of the map, they will be ineffective against the mix- 
ing of power due to a point source mask. We propose here 
a method for dealing with such issues, which in essence 
is the following: 



A. Prewhitening of noiseless CMB maps 

Let us denote the pure CMB temperature map as T{9). 
An important feature of the power spectrum of the CMB 
C( is that beyond a multipole of ~ 1000, it is approxi- 
mately proportional to £~ A . Such sharp fall in the power 
makes it extremely prone to aliasing of power across 
multipoles due to a point source mask. If we perform 
an operation akin to taking the Laplacian of this map, 
then we would effectively multiply the power spectrum 
by £ A on all scales, thereby making the processed power 
spectrum nearly white over the large £ tail. This would 
thereby minimize aliasing of power. In the following we 
propose a method of achieving this, by a combination 
of two operations which we call "disc-differencing" and 
"self-injection" . 



1. Disc-differencing 

If we convolve the map with a circular disc of radius 
R, generating, 



t r {6) = J d 2 e'w R {o' - e)T{e') 

where Wr is the top-hat filter, 



W R {6) 



otherwise, 



(40) 



(41) 
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FIG. 6: Comparison of the probability distributions of the estimated power spectrum at three points (pixels) in the two- 
dimensional Fourier (I) space. From left to right, these points are (fj 2 ) = (630,630), (1980,630) and (2880,630), for which 
the modulus £ has been indicated in each figure. Upper Panel: Power spectrum estimation with N res = 2.0 and Nt ap = 3 2 . 
The open circles (black), the diamonds (red) and the triangles (green) indicate the probability distribution of the quantity 
£(£ + l)C e /(2n) as estimated via the AMTM, MTM and PM methods from 5000 Monte Carlo simulations, respectively. The 
respective approximate theoretical forms as discussed in § HE are over-plotted as the continuous curve (black), the dotted 
curve (red) and the dashed curve (green) for each of the methods. The mean degree of freedom of the chi-square for each 
method is also indicated as va for AMTM, vm for MTM and vp for the periodogram, PM. Each curve is also accompanied by 
a vertical line of the same style (and color) representing the mean value obtained from the Monte Carlo simulations. In each 
figure, the continuous (black) vertical line corresponding to the mean of the AMTM method, is also the value closest to the 
true power spectrum. It is actually the unbiased value of the pseudo power spectrum (see § V). Lower Panel: Same as above 



but with iV r 



3.0 and Ntap = 5^ 



then in Fourier space, we will have 

Tr(£) = W R (£)T(£) (42) 
where the Fourier space window, W(£) is given by, 

MR) 



W R (£) = 2- 



IR 



(43) 



Now let us consider smoothing the map with another 
top-hat window of radius 3i?, giving, 



t 3R (0) = J d 2 e'w 3B xe' -e)T(e'). 

We then take the difference map, 

Trsr = T R - T 3R , 

which in Fourier space reads, 

T R - 3R (£) = W R - 3R {£)T(t) 

'Ji{lR) ,h{UR) 



(44) 



(45) 



= 2 



IR 



UR 



(46) 

T{1). (47) 



Remembering the asymptotic expansion of J\{x) for 
x « a/2, 



i 



(48) 



which shows that for values of £ such that £R < a/2, we 
have, 

T R - 3R (l) ~ (£R) 2 T(£), (49) 

which implies that the power spectrum is now, 

(2TrfC*- 3R ~ (T R _ 3R (£)*T R _ 3R (£)) = (2tt) 2 £ 4 R 4 C £ , 

(50) 

which is the desired form. The disc-difference window, 
W R ^ 3R (£R) is shown in Fig. 9. Beyond I = y/2/R, the 
function starts falling again. Thus by choosing the radius 
R judiciously, one can prewhiten a desired range of the 
power spectrum. 

An important aspect of this method is that it is an ef- 
fective C" 1 / 2 operation on the map (C is the covariance 
matrix), which is manifestly local, and therefore does not 
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FIG. 7: Comparison of the probability distributions of the estimated power spectrum in bins. Each panel is for a different 
combination of the JVres and Nt ap as indicated on the top of the middle figure. The open circles (black), the diamonds (red) 
and the triangles (green) indicate the probability distribution of the quantity lit + 1)Ci/(2tt) in each bin as estimated via 
the AMTM, MTM and PM methods from 5000 Monte Carlo simulations, respectively. The respective approximate theoretical 
forms as discussed in § HE (see eq. (37)) are over-plotted as the continuous curve (black), the dotted curve (red) and the 
dashed curve (green) for each of the methods. The mean degree of freedom of the chi-square for each method is also indicated 
as va for AMTM, vm for MTM and vp for the periodogram, PM. The number of pixels in each bin is also indicated as Nb. 
Note the the periodogram is absent in each of the rightmost plots, as it is highly biased and lies outside the range plotted. 



suffer from effects due to edges, which is a common prob- 
lem in case of the full Fourier space operation. Maximum 
likelihood methods of estimation of power spectra nat- 
urally involve the C _1 operation [30, 31]. However, for 
high resolution experiments, the numerical costs for com- 
puting the C _1 matrix can be prohibitively large. Also, 
such an operation is non-local and mixes modes from 
masks with the map in a non-trivial way. The method 
proposed here is approximate but is simple to implement 
as it involves only convolutions in real space: as such, 
its effect can be quantified, propagated or undone very 
easily. 

Figure 10 illustrates the effect of the disc-differencing 
operation on the CMB power spectrum, for a radius 



R = V . This means that the prewhitening window turns 
around at I ~ V~2/R — 5000. The application of the 
disc-differencing produces the processed power spectrum 
shown by the dashed curve in the figure, which has much 
smaller dynamic range than the original one shown by 
the dotted curve. 



2. Self-injection 

One undesirable effect of the application of the disc- 
differencing function is that it makes the lower multipolc 
part of the CMB power spectrum {I < 1000) a steeply 
rising function, which may aggravate aliasing of power. 



13 



0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 0.0 2 0.4 0.6 0.8 1.0 

fj = 2.0; N a = 3 N = 3.0; N a = 5 N = 4.0; N a = 4 




1000 2000 3000 4000 5000 6000 1000 2000 3000 4000 3000 6000 1000 2000 3000 4000 5000 6000 



FIG. 8: Covariance matrix of the bandpowers estimated via AMTM for three different parameter settings; from left to right 
these are: N re3 = 2.0, Nt ap = 3; N res = 3.0, Nt ap — 5 and N re3 = 4.0, Nt ap = 4. Each square in the image represents a 
bin in I. For the N res = 2.0 and JV res = 3.0 cases, we chose the bin-widths to be twice the fundamental resolution element, 
i.e. A£ — 21 f un d = 180, while for the N res = 4.0 case it was taken to be 4£f un d = 360. There is appreciable covariance only 
between bins inside the resolution A£w = 2 N rea £f un d set by the taper, and the covariance drops drastically beyond that 
frequency. 
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FIG. 9: "Disc-difference" function Wrsr discussed in the 
text (solid line). The dashed line represents the function x 4 . 



A simple way to deal with this problem, is to add a small 
fraction of the original map back into the disc-differences 
map, a process we call "self-injection" . If a fraction a <C 
1 of the map is self-injected after the disc-differencing, 
the processed power spectrum looks flat at all multipoles, 
and is conveniently given by a multiplication of the true 



power spectrum with an analytic function, 

C[ w = (W R . 3R (m) + a) 2 C e . (51) 

This is illustrated for a = 0.02 in Fig. 10. 

Having laid out the basic theory of prewhitcning, we 
will now describe a concrete example of power spectrum 
estimation of a CMB map with a point source mask to 
judge the effectiveness of this method in recovering the 
true power spectrum. To this end we simulate a 8° x 8° 
map as a Gaussian random realization from a theory 
power spectrum. The map has 768 pixels on a side, mak- 
ing the pixel scale A9 = 0.625'. We cut out a 4° x 4° 
section from this map to impose non-periodic boundary 
conditions. We then simulate a point source mask for 
this smaller map, which is unity everywhere except in- 
side randomly positioned discs of various radii, where it 
is zero. 

To motivate the necessity of prewhitening, we first ap- 
ply the mask directly to the map and take its AMTM 
power spectrum with N res = 3.0 and N tap = 3 2 . The re- 
sult is shown in Fig 11. Note that we expressly use only 
the first few tapers with least leakage to ensure that mode 
coupling is minimized due to sharp edges of the map. 
However, this choice has no bearing against the aliasing 
of power due to holes in the map and as expected, we find 
a lot of power aliased from the low to the high multipoles. 
Next we perform a disc-difference operation on the orig- 
inal map with R = 1', followed by a self-injection with 
a = 0.02. Then, we apply the mask on this prewhitened 
map and perform AMTM power spectrum estimation on 
it. By comparing with the expected theoretical power 
spectrum, wc find that aliasing is almost non-existent in 
the power spectrum of the prewhitened and masked map. 
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FIG. 10: Effect of disc-differencing and self-injection on the 
power spectrum. The dotted line is the true power spectrum 
with a large dynamic range. The disc-differencing operation 
alone produces the dashed curve which has a much smaller dy- 
namic range, but is steeply rising at low multipoles. The disc 
radius used is R — 1'. The dot-dashed curve shows the true 
power spectrum multiplied by a constant a 2 where a = 0.02. 
If we disc-difference the map followed by self-injection of a 
fraction a of the map, then the power spectrum of the pro- 
cessed map is the solid curve (given by (51)) which is conve- 
niently flat over the range of multipoles. 
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FIG. 11: Prewhitening and AMTM as a remedy to aliasing 
of power due to point source mask. The dotted curve repre- 
sents the input power spectrum from which the map is gen- 
erated. The triangles represent the recovered AMTM power 
spectrum (N res — 3.0, Ntav ~ 3 2 ) of the map after a point 
source mask is applied directly to it. If, on the other hand, 
the map is first prewhitened (see text) and the mask is sub- 
sequently applied, one obtains the diamonds as the AMTM 
power spectrum. The solid curve is the theoretical predic- 
tion for the prewhitened power spectrum. Thereafter, one 
divides the spectrum by the prewhitening transfer function 
(see (51), obtaining a nearly unbiased estimate, denoted by 
the open circles. Note that the AMTM power spectra appear 
smoother than the true spectra as the former is convolved 
with the window function of the taper. 



Then we simply divide this power spectrum by the an- 
alytical transfer function for the prewhitening operation 
(51), to obtain a nearly unbiased estimate of the true 
power spectrum. 



B. Prewhitening of Noisy Maps 

In a real experiment, the map will be convolved with 
the instrument beam and will contain noise from the in- 
strument as well as other astrophysical signals. To simu- 
late the simplest of such situations, we convolve the map 
from the previous step with a Gaussian beam of full- 
width- at-half-maximum (FWHM) of 5'. Next, we add 
Gaussian white noise at a level of 3 /iK per sky pixel 
(defined as an area of FWHM 2 on the sky) . The power 
spectrum of the map so generated has the usual initial 
sharp fall with multipole £, but then flattens out beyond 
the multipole where white noise starts to dominate. One 
immediate consequence of this is that the disc differenc- 
ing operation that multiplies the power spectrum by £ A 
will make the white noise part rise with £ instead of 



remaining flat. One can try to minimize this effect by 
judiciously choosing the disc radius R so that the disc- 
differencing window function turns over at the value of 
£ where white noise starts to dominate. However, there 
may be components to the power spectrum other than 
white noise which may rise faster than the fall of the 
disc-differencing window. All these cases can be effec- 
tively dealt with by introducing another component to 
the prewhitening operation: namely a convolution with a 
Gaussian window after the disc-differencing and the self- 
injection steps. If we convolve the map with a Gaussian 
of FWHM 6 S then after evaluating the power spectrum 
we can remove its effect by dividing it by the multipole 
space form cxp(£(£ + l)0 s /(8 log 2)). In Fig. 12 we illus- 
trate the prewhitening in presence of white noise. 



V. MODE-MODE COUPLING AND 
DECONVOLUTION 

As discussed earlier, an undesirable feature of the Mul- 
titaper method is the loss of spectral resolution. For any 
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FIG. 12: Same as in Fig. 11, but for a map with white noise. 
In addition to the standard prewhitening operation, a Gaus- 
sian smoothing has been applied to flatten the tail of the 
prewhitened map. 



N res > 1, a Multitapcr estimator smooths the power 
spectrum with a frequency-space window which is wider 
than the fundamental resolution set by the size of the 
map. For power spectra which are highly structured, this 
poses the problem of diluting interesting features which 
may render the power spectrum less useful as a direct 
probe of the underlying phenomenon. For example, in 
case of the CMB, (A)MTM leads to the smoothing of 
acoustic features. 

In CMB analyses, the problem of recovering the true 
power spectrum from one that has been convolved with 
a window function is a well studied problem. It arises 
naturally in full-sky CMB experiments because of the 
application of point source and galactic masks. A de- 
tailed account of the procedure to deal with mode-mode 
coupling in the spherical harmonic space can be found 
in [15]. In what follows, we will discuss the mode-mode 
coupling in flat space and develop the method of decon- 
volving the effective tapering window function from the 
power spectrum. 

Consider a homogeneous, isotropic and zero-mean 
Gaussian random process realized on a map, T(x), with 
a power spectrum P(k), 



(T(k)*T(k')) = (27r) 2 (5(k- k')P(fc). 



(52) 



If the map is multiplied by a window W(x), then the 
Fourier transform of the resulting map is given by the 
convolution, 



T w (k) 



d 2 k' 

(2tt) s 



T(k')W(k-k'). 



(53) 



In the ensemble average sense, we therefore have, 
1 



iV(k) 



(2tt)s 



(T w *{k)T w (k)) 



(54) 



d 2 k' \W(k-k')\ 2 
(2^)2 (2wf 



P(k'), (55) 



which is often referred to as the pseudo power spectrum 
in CMB literature. All the multitapered power spec- 
tra discussed so far are pseudo-power spectra in this 
sense, and comparing with (30), it is easily seen that 
for AMTM, the power spectrum of effective window that 
couples to the true power spectrum is given by, 



W 



AMTM 



(k) 



V b 2 (k)V( a 



ia 2 ;JV 1 JV 2 ) 



(k) 



T bf Jk) 



(56) 

One is obviously interested in the angle averaged pseudo 
power spectrum, 



P w (k) = 



d6 
W) 



P w {k) 



(57) 



and as demonstrated in appendix B, this can be expressed 
in terms of the true power spectrum as, 



P w (k)) = / dk'M(k,k')P(k') 



(58) 



where M is the mode-mode coupling kernel and depends 
on the power spectrum of the window W. After the 
power spectrum is binned, the above relation can be con- 
veniently expressed as a matrix multiplication (see ap- 
pendix B,), 



P 



(59) 



where b are the indices for the bins and M is the binned 
mode -mode coupling matrix. The binned mode-mode 
coupling matrix is more stable to inversion than the un- 
binned one which tends to be nearly singular, and there- 
fore an estimate of the true power spectrum can be ob- 
tained by inverting the above relation, 



P b = (MT 1 ) W P™ . 



(60) 



We first illustrate the deconvolution of the power spec- 
trum without a point source mask. For this purpose, 
we simulate a large 16 degree square map, with 768 pix- 
els on a side, using a theoretical power spectrum. We 
then perform power spectrum analysis of the central 8 
degree square portion of it. For the periodogram, wc 
simply zero out the portion of the larger map outside 
the central 8 degree square, so that the mask becomes 
a zero-padded top-hat. For the AMTM, we discard the 
outer regions and perform the adaptive multitaper anal- 
ysis with N res = 3, N ta p — 5 2 on the central 8 degree 
square. We generate the mode coupled power spectrum 
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FIG. 13: Deconvolution of the power spectrum. Left panel: Deconvolution of the periodogram with top-hat weighting. The 
black squares represent the binned power spectrum obtained directly from the map using the periodogram (straight FFT) 
method. As discussed in the text, this power spectrum is the true power spectrum convolved with the mode-mode coupling 
matrix due to the top-hat, the theoretical expectation for which is displayed as the blue line. The red points represent the 
binned power spectrum deconvolved via (60). The black line is the input power spectrum. The deconvolved binned power 
spectrum is to be compared with the binned input power spectrum which is displayed as the black histogram. All points 
displayed are the mean of 800 Monte Carlo simulations and the error bars correspond to the 2a spread in their values. Right 
Panel: Same as above, but for the AMTM method. The mode-coupled power spectrum and the corresponding theoretical 
curve in this case have been artificially shifted below the deconvolved power spectra for easy viewing. As discussed in the text, 
the mode-coupled power spectrum produced by the AMTM method is a nearly unbiased estimate of the true power spectrum, 
while the mode-coupled periodogram (left panel) is highly biased at large multipoles. This bias causes the error bars in the 
deconvolved periodogram to be much larger than those in the deconvolved AMTM power spectrum at large £, as shown in 
Fig. 14. 



with each method and then deconvolve it using (60). The 
results from 800 Monte Carlo experiments are shown in 
Fig. 13. 

Firstly, we note that the construction of the mode- 
mode coupling matrix, as delineated in appendix B, en- 
ables us to accurately predict the mode-mode coupled 
theoretical power spectrum. As such, likelihood analyses 
for cosmological parameters can potentially be performed 
with the mode-coupled AMTM power spectrum, and this 
is specially appealing because the statistical properties of 
the mode-coupled binned power spectra can be precisely 
predicted (see § HE). 

Secondly, we find that the deconvolution of the peri- 
odogram as well as the AMTM power spectrum produces 
an unbiased estimate of the input power spectrum, but 
with a very important difference: the uncertainty in the 
deconvolved periodogram at £ > 1500 is larger than that 
in the deconvolved AMTM spectrum, and the difference 
continues to grow to factors of several for higher mul- 
tipoles. This is seen clearly in Fig. 14, where we have 
plotted the fractional errors in the deconvolved binned 
power spectrum against the bin centers. For instance, at 
I ~ 3000 the periodogram method produces ~ 3 times 
larger error bars. This owes to the fact that due to spec- 



tral leakage, the periodogram produces a pseudo power 
spectrum that is highly biased relative to the input the- 
ory over these higher order multipoles. This bias adds 
to the uncertainty in the deconvolved power spectrum. 
On the other hand, AMTM produces a nearly unbiased, 
albeit smoothed, pseudo power spectrum which, when 
deconvolved, does not incur any excess uncertainties. Ig- 
noring point source masks for the moment, which are 
very specific to the CMB, this result alone immediately 
elevates the multitaper method to a far superior status 
than the periodogram, as a power spectrum estimation 
method for highly colored spectrum measured from finite 
maps. 

One must bear in mind, that over the range of mul- 
tipoles where leakage is not a serious problem, and the 
periodogram is essentially unbiased, the errors bars ob- 
tained from the periodogram are the smallest. This is 
because the periodogram makes use of all the modes that 
are available from the entire map, while each taper in a 
MTM process makes use of a certain fraction of the map. 
By using several tapers, one can compete with the pe- 
riodogram error bars over the nearly white part of the 
spectrum. This is essentially seen in the convergence of 
the fractional error bars from the two methods in the low 
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FIG. 14: Fractional errors in the deconvolved binned power 
spectrum. The open circles represent the fractional errors for 
the periodogram method (left panel of Fig. 13). The filled cir- 
cles represent the same for the AMTM method ( right panel 
of Fig. 13). Although the deconvolved power spectrum ob- 
tained from either method is an unbiased estimate of the true 
power spectrum, the errors from the periodogram method are 
much larger at high £ because of the highly biased nature of 
the mode-coupled periodogram at those multipoles. 



£ portion of Fig. 14. This, however, is not in contradiction 
to the fact, as shown in § HE that it is possible to have 
lower uncertainties in any bin for the pseudo power spec- 
trum with AMTM than with the periodogram, by choos- 
ing a high enough N res (and consequently high N tap ). 
Unfortunately, this advantage goes away when the power 
spectrum is deconvolved. 

Next we turn to the practical issue of dealing with 
a point source mask for a CMB map. As already dis- 
cussed in § IV, the holes in the mask couple power 
over various scales, and neither the periodogram nor 
the AMTM remain an unbiased method at large mul- 
tipoles, necessitating the prewhitening of the map. To 
study the effect of the mask and prewhitening on the de- 
convolved power spectrum, we multiply our map with a 
point source mask and repeat the Monte Carlo exercise 
with and without prewhitening. The mask has a cover- 
ing fraction of 98% and the holes are of random sizes. 
For the prewhitening operation we use the parameters 
R = V and a = 0.02. The results on the fractional er- 
rors for this case are displayed in Fig. 15. This figure 
shows the power of prewhitening approach. Prewhiten- 
ing significantly reduces the errors. It is obvious that 
like the periodogram, AMTM is defenseless against the 
point source holes in the mask. This is expected because 
a taper with holes no longer retains the nice property of 
being a band-limited function. For both the top-hat and 
the effective tapered window, the sidclobes become dom- 



FIG. 15: Effect of prewhitening on the errors in the decon- 
volved power spectrum in presence of a point source mask. 
Top Panel: For the AMTM method: The open circles rep- 
resent the fractional la error bars on the deconvolved power 
spectrum when prewhitening is not performed, showing that 
holes in the point source mask render the AMTM method 
biased at high multipoles and lead to large error bars. The 
filled circles represent the same after prewhitening has been 
performed, showing that prewhitening remedies the leakage 
of power and makes the power spectrum estimator nearly 
unbiased. Bottom Panel: The same as above, but for the 
periodogram or straight FFT. Note that prewhitening, if per- 
formed properly, makes the periodogram as good a power 
spectrum tool as the AMTM. 



inated by the power produced by the holes. In such a 
situation, prewhitening the map becomes a necessity to 
reduce the uncertainties in the deconvolved power spec- 
trum. With proper prewhitening both the periodogram 
and the AMTM method can produce nearly unbiased 
pseudo power spectra and consequently small error bars 
at the high multipoles, as evident from Fig. 15. It is im- 
portant to note here that the power spectrum of a top- 
hat has greater sidelobc power than the effective taper 
in the intermediate range of multipoles between the re- 
gions where the central lobe is falling and where the point 
source power starts to dominate. Therefore, depending 
on the quality of prewhitening that can be performed on 
a map, AMTM may be a more reliable method to apply 
on the prewhitencd and masked map, rather than the 
periodogram. 



VI. CONCLUSION 

Power spectrum estimation on small sections of the 
CMB sky is a non-trivial problem due to spectral leak- 
age from the finite nature of the patch, which is further 
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compounded by the application of point source masks. 
The direct application of standard decorrelation tech- 
niques, like the MASTER algorithm [15], to obtain an 
unbiased estimate of the power spectrum leads to unnec- 
essarily large uncertainties at high multipoles due to the 
highly biased nature of the pseudo power spectrum at 
those multipoles. We have put forward two techniques 
to reduce the uncertainties in the deconvolved power 
spectrum. First, we have formulated a two-dimensional 
adaptive multitaper method (AMTM) which produces 
nearly unbiased pseudo power spectra for maps without 
point source masks, by minimizing the leakage of power 
due to the finite size of the patch. This is achieved at 
the cost of lowered spectral resolution. The deconvolu- 
tion of the pseudo power spectrum so produced, leads 
to an unbiased estimate of the true power spectrum 
that has several times smaller error bars at high multi- 
poles than the deconvolved periodogram. In presence of 
point source masks, however, this method becomes non- 
optimal because the pseudo power spectrum estimated 
even by AMTM is no longer unbiased. To deal with the 
point source mask, we have put forward a novel way of 
prcwhitening a CMB map, with manifestly local oper- 



ations which has simple representations in the Fourier 
space. This operation produces a map, the power spec- 
trum of which has several orders of magnitude lower dy- 
namic range than the original map. This renders the 
leakage of power due to holes and edges a relatively be- 
nign issue for the prewhitcned map. If the prcwhitening 
operation can be tuned to make the power spectrum of 
the map nearly white, a pseudo power spectrum obtained 
via a simple periodogram may be nearly unbiased and 
therefore, can be deconvolved to give a precisely unbiased 
estimate of the power spectrum, thereby avoiding unnec- 
essarily large error bars at large multipoles. If the map 
cannot be made sufficiently white for a periodogram, an 
AMTM method can be applied to the prewhitened map 
to guard against leakage and achieve the same result. We 
have shown that by applying these methods, one can re- 
duce the error bar in the small scale power spectrum by 
a factor of ~ 4 at i ~ 3500. This can be translated into 
a many-fold reduction in the required integration time of 
a CMB experiment to achieve some target uncertainty 
on the small scale power spectrum than that dictated by 
standard techniques. 



I 

APPENDIX A: STATISTICAL PROPERTIES OF MULTITAPER ESTIMATORS 



The statistical properties of the multitaper spectral estimators follow from the basic result (32) for the distribution 
of the periodogram. In the following, we will use the equivalent degrees of freedom argument to derive the distribution 
functions for the MTM estimators. The first result we will use is that the probability distribution of an eigenspectrum 
p( QlQ2 )(k) at each frequency-space point k (i.e. at each pixel) will have the same form as (32) for asymptotically 
large N. Since the final spectral estimate is a weighted sum of these eigenspectra, it is reasonable to assume that the 
former is distributed also as a scaled \ 2 variable. Let us hypothesize, 

pAMTM^d^ (A1) 

where both v and a are unknown. Now, we will make use of the facts that ^p AM ™^ = (ax 2 ,) = av and var(P AMTM ) = 
var(ax^) = 2a 2 v. Therefore, 

9^pAMTM^ 2 /pAMTM\ 
V = var(P AMTM ) 5 a= v ' (A2) 

Since the eigenspectra P"i Q 2 arc asymptotically uncorrclated and are also unbiased estimators of P(k), we obtain, 
using the definition of the AMTM estimator (30), 

(P AMTM ) = P(k), (A3) 
var(P AMTM ) = P 2 (k) ^ ai ' a2 (ai " 2)V 2 (A4) 



V b 2 fkl 



for N — > oo, which immediately gives us 



, ( k) = ^ ni "; ttia y . (as) 



which is equation (34). 
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Now we turn to the approximate form of the distribution for the binned power spectrum, P B (35). We again 
postulate that P B is distributed as cx„ so that, 

2(P*) 2 (PB) 

^-^n - c -—- (A6) 

Now note that because the images or maps we will be concerned with are real, only the half-plane in the k space are 
independent because of reflection symmetry i.e P(k) = P(— k). For the time being, let us assume that all k space 
pixels are uncorrelated, so that all pixels in one half of the k plane are independent of each other. Later, we will 
correct this for the fact the tapers highly correlate adjacent pixels in each resolution element. With these assumptions, 
(35) can be re-written as, 

pB ^ = W £ P AM ™(h,ki). (AT) 

b i,jeb;kj>0 

Therefore, the variance is given by, 

™(p B (h)) = ± £ ^(^ AM ™ft-y) 2 



4 / nAMTM/i' ' . » 2 



(p AMTM (|k| = h)y 



i,jeb 

where in the last step we have assumed that the power spectrum is slowly varying inside a bin. The factor of 2 inside 
the summation goes away because we have extended the summation to include the lower half plane. Using (A6), we 
therefore get 

1 2 sr- 1 

Now, we will correct for the fact that all pixels in the upper half plane are not independent. In fact for a taper with 
a given N res , approximately N^ es pixels get highly correlated. Therefore we can think of "super-pixels" of dimension 
N res x N res , which arc independent of each other. If we repeat the calculation above with this assumption, then we 
obtain the result, 

-j res \ 



(A10) 

i,jeb 



which essentially reduces the degrees-of-freedom by the number of pixels that fall in each super-pixel. This is equation 
(37). 

APPENDIX B: MODE-MODE COUPLING MATRIX 

In the following, we derive the form of the mode-mode coupling kernel. For the sake of completeness, we repeat 
and expand upon some of the calculations from appendix A.l of [15]. Note that the symbol k (frequency) in the said 
appendix of that paper is smaller by a factor of 2-7T than the k (angular frequency) used in the present section. In 
other words, the k used here is readily identified with the multipole £ in the flat-sky approximation, while it would 
stand for £/(2n) in the aforementioned appendix. 

We start with the definition (57) of the angle averaged pseudo power spectrum and with the aid of (54) express it 
as, 

We can express the Fourier mode of the window appearing above as 



W(kx - k a ) = / d z k 3 W(k 3 )6 z (k 3 - ki + k 2 ) (B2) 
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which immediately gives, 



|W(ki - k 2 



2 



2tt J dk 3 k 3 W{k 3 )6 2 (k 3 - kx + k 2 ) (B3) 



(2tt)2 

where we have introduced the power spectrum of the window. 

(2nfW(k) = J -^L W *(k)W(k). (B4) 

Substituting in (Bl), we get, 

[ pw ^) = J^sJ / dd2P[h2) I dh w^ 2 ^ 3 ki + k2 ) (B5) 

dk 2 M klk2 P{k 2 ) 7 (B6) 



M klk2 = 77r~r / dk 3 k 3 W{k 3 )J{k u k 2 ,k 3 ) (B7) 



where we have introduced the mode-mode coupling matrix, 

fc 2 

(27T) 



with, 



(27T) 

It is instructive to compare (B7) with the full-sky result from appendix A. 2 of [15] 



J(h, k 2 ,k 3 ) = 1-^-1 d0 2 S 2 (k 3 - ^ + k 2 ) (B8) 



Mi 



2£ 2 ' 



1<2 



^rE^ + ^f^ o 2 o ) • ( B9 ) 



It can be shown that for large 



*i 4 la x 




J(£ul 2 ,£ 3 ) (BIO) 



so that if we identify k{ with the large £ limit of the spherical mode-mode coupling matrix goes correctly to the 
flat-sky expression (B7). Note that the appearance of an extra factor of 2 in the large £ limit of the summand in (B9) 
over the integrand in (B7) is fine because the sum is restricted to values £ 3 which make {£\ + £ 2 + £ 3 ) even, while the 
integral over k spans the entire allowed range for each (ki, k 2 ) pair. 

We now turn to the evaluation of the J function. Using the factorized form of the delta function in plane polar 
coordinates, <5(r — r') = S(r — r')8{6 — 6')/r , we can write, 



,k 2 ,k 3 ) = [d6 2 6 2 (ki - |k3+k 2 |)/*i. (Bll) 
(A) J 



The integral over 9 2 can be performed by using the following property of the delta function, 

'taCM-ElFfcyf (B12) 

where Xi are the roots of g(x) = 0. In our case there are two roots and the integral finally yields 

■7(fci,fc 2 ,fc 3 ) = - ; 1 (B13) 

7T y/-K{K - 2ki){K - 2k 2 )(K - 2k 3 ) 

for \k 2 — k 3 \ < ki < k\ + k 2 , and zero outside the interval, where K = ki + k 2 + k 3 . As a word of caution against 
using this form directly into the integral (B7), we would like to point out that due to the diverging nature of the 
integrand near the edges of the allowed interval, numerical integration schemes that span the entire range will be 
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unreliable. Instead wc find, drawing parallels with the spherical harmonic case, that approximating the integration 
in (B7) by a sum over integral values of fc, by interpolating W(fc) onto integers, provides a very stable and reliable 
way of computing the mode coupling matrix. In fact, the boundaries of interval can be included in the sum by using 
the large t limit of the 3 — j functions at those points, 



J y_^L for £ 3 = 



(B14) 



At the end of the exercise, we therefore generate the mode-coupling matrix at integer subscripts, which by abusing 
notation a bit, we will denote by Mw. 

Next, we consider the binning of the pseudo power spectrum. In order to construct the mode-mode coupling kernel 
relevant to a binned power spectrum, we first define the binning operator B, which performs the binning of the integral 
indexed quantities to binned values and its reciprocal operation U. Given the lower and upper boundaries, and 
^high °f a b m b, a simple form of these operators can be written as, 

{ 1° if o < P b < P <r P b 

W^U> *2<« IOW <*<* high (B15) 
0, otherwise, 

and 

{— if 9 < P b < P <r f b 
js, a^S^SK^ (B16) 
0, otherwise. 

Here, a is chosen to make the power spectrum "flatter" and for the damping tail of the CMB, a suitable value is 
a = 4. 

Note that the pseudo power spectrum realized on the two-dimensional k-space may be directly binned into the bins 
b, by averaging the value in the pixels that fall inside the annuli demarcated by the bin boundaries, without having 
to go through an intermediate step of interpolating the power spectrum onto integers. The binned pseudo power 
spectrum is, therefore related to the binned true power spectrum as 

(PF) = BuMwPf 

= BbiMu>Ue>b> Bbi£iPt' 

= M w Pb>, (B17) 
where we have defined the binned mode-mode coupling matrix as M = IE 



ACKNOWLEDGMENTS 

We dedicate this paper to the memory of eminent geo- 
physicist F. Anthony Dahlen (1942-2007), who had sug- 
gested to us the idea of multitapering for the CMB. We 
would like to thank Frederick J. Simons, Mark A. Wiec- 
zorek, Neelima Sehgal and Tobias A. Marriage for en- 
lightening discussions and Viviana Acquaviva for useful 
comments on the manuscript. SD would like to acknowl- 



edge the warm hospitality extended by the mentors and 
colleagues at Jadwin hall, specially Lyman Page and Jo 
Dunklcy, while the Department of Astrophysics was un- 
der renovation. AH would like to thank Eric Hivon for 
useful discussions. AH is supported by the LTSA pro- 
gram. SD is supported by the Charlotte Elizabeth Proc- 
ter Honorific Fellowship from Princeton University and 
NSF grant 0707731. DNS acknowledges support from 
NASA ATP grant NNX 08AH30G. 



[1] A. Hajian and T. Souradeep, Phys. Rev. D 74, 123521 

(2006), arXiv:astro-ph/0607153. 
[2] E. Komatsu, A. Kogut, M. R. Nolta, C. L. Bennett, 

M. Halpern, G. Hinshaw, N. Jarosik, M. Limon, S. S. 

Meyer, L. Page, et al., Apjs 148, 119 (2003), arXiv:astro- 

ph/0302223. 

[3] K. M. Gorski, A. J. Banday, C. L. Bennett, G. Hinshaw, 



A. Kogut, G. F. Smoot, and E. L. Wright, ApJ 464, 
L11+ (1996), arXiv:astro-ph/9601063. 

[4] M. R. Nolta, J. Dunkley, R. S. Hill, G. Hinshaw, E. Ko- 
matsu, D. Larson, L. Page, D. N. Spergel, C. L. Bennett, 

B. Gold, et al., ArXiv e-prints 803 (2008), 0803.0593. 
[5] G. Hinshaw, J. L. Weiland, R. S. Hill, N. Odegard, 

D. Larson, C. L. Bennett, J. Dunkley, B. Gold, M. R. 



22 



Greason, N. Jarosik, et al., ArXiv e-prints 803 (2008), [20 
0803.0732. 

[6] C. L. Reichardt et al. (2008), 0801.1491. 

[7] QUaD collaboration: C. Pryke, P. Ade, J. Bock, M. Bow- [21 
den, M. L. Brown, G. Cahill, P. G. Castro, S. Church, 
T. Culverhouse, R. Friedman, et al., ArXiv e-prints 805 [22 
(2008), 0805.1944. 

[8] http : //www. physics .princeton. edu/ act/. 

[9] http://pole.uchicago.edu. 
[10] http://www.rssd.esa.int/Planck/. [23 
[11] D. B. Percival and A. T. Walden, Spectral Analysis for 
Physical Applications, Multitaper and Conventional Uni- 
variate Techniques (Cambridge Univ. Press, New York, [24 
1993). 

[12] M. Tegmark, Phys. Rev. D 56, 4514 (1997). [25 
[13] O. Dore, R. Teyssier, F. R. Bouchet, D. Vibert, 

and S. Prunet, A&A 374, 358 (2001), arXiv:astro- [26 

ph/0101112. 

[14] K. M. Smith, O. Zahn, and O. Dore, Phys. Rev. D 76, [27 
043510 (2007), arXiv:0705.3980. 

[15] E. Hivon, K. M. Gorski, C. B. Netterfield, B. P. Crill, [28 
S. Prunet, and F. Hansen, Astrophys. J. 567, 2 (2002), 
arXiv:astro-ph/0105302. [29 

[16] P. Welch, Audio and Electroacoustics, IEEE Transac- 
tions on 15, 70 (1967), ISSN 0018-9278. [30 

[17] D. Thomson, Proceedings of the IEEE 70, 1055 (1982). 

[18] D. Slepian, AT T Technical Journal 57, 1371 (1978). [31 

[19] B. Bell, D. B. Percival, and A. T. Walden, Journal of 

Computational and Graphical Statistics 2, 119 (1993), [32 
ISSN 10618600, URL http://www.jstor.org/stable/ 
1390958. 



H. Cramer, The Annals of Mathematics 41, 215 (1940), 
ISSN 0003486X, URL http://www.jstor.org/stable/ 
1968827. 

J. L. Doob, SIAM Review 5, 172 (1963), URL http: 
//link . aip . org/link/?SIR/5/172/2. 
M. B. Priestly, Journal of the Royal Statistical Soci- 
ety. Series A (Statistics in Society) 151, 573 (1988), 
ISSN 09641998, URL http://www.jstor.org/stable/ 
2983035. 

Barr, D. L., Brown, W. L., and Thompson, D. J., Le 
Journal de Physique Colloques 49, C6 (1988), URL 
http : //dx . doi . org/doi/10 . 1051/ jphyscol : 1988630. 
T. Bronez, Acoustics, Speech and Signal Processing, 
IEEE Transactions on 36, 1862 (1988), ISSN 0096-3518. 
T.-C. Liu and B. van Veen, Signal Processing, IEEE 
Transactions on 40, 578 (1992), ISSN 1053-587X. 
A. Hanssen, Signal Process. 58, 327 (1997), ISSN 0165- 
1684. 

M. A. Wieczorek and F. J. Simons, Geophysical Journal 
International 162, 655 (2005). 

F. A. Dahlen and F. J. Simons, ArXiv e-prints 705 
(2007), 0705.3083. 

F. J. Simons, F. A. Dahlen, and M. A. Wieczorek, ArXiv 

Mathematics e-prints (2004), math/0408424. 

M. Tegmark, Phys. Rev. D 55, 5895 (1997), arXiv:astro- 

ph/9611174. 

S. P. Oh, D. N. Spergel, and G. Hinshaw, Astrophys. J. 
510, 551 (1999), arXiv:astro-ph/9805339. 
sine refers to the sinus cardinus i.e. sinx/x. Power spec- 
trum of a two-dimensional top-hat window is a product 
of two sine 2 functions. 



